Skip to main content

Programación y Análisis de Datos Geoespaciales

1 bibliography: [book.bib, packages.bib]

Placeholder

1.1 DESCRIPCIÓN DEL CURSO

1.1.1 OBJETIVOS

1.2 Evaluación

1.2.1 tareas

1.2.2 proyectos

1.2.3 calificaciones

1.3 Calendario

1.3.1 Módulo 1: Introducción a R y RStudio

1.3.2 Módulo 2: R Markdown & Github

1.3.3 Módulo 3: Manipulación de Datos con Tidyverse

1.3.4 Modulo 5: Visualización de Datos

1.3.5 Módulo 6: Introducción al Análisis Geoespacial

1.4 Recursos adicionales

2 Introducción

Placeholder

2.1 Objetos

2.2 Variables

2.2.1 Tipos de variables

2.3 Funciones

2.4 Vectores

2.5 Instalar librerías

2.6 R Notebook

2.7 Leer datos

2.8 Ejercicios

3 Principios de Tidydata

Placeholder

3.1 Definiciones

3.2 Lógica y funciones

3.2.1 Funciones esenciales

3.2.2 pipe

3.2.3 Agrupar

3.2.4 Filtrar

3.2.5 Crear nuevas variables

3.2.6 Responder:

3.3 Cambiar formato de tabla

3.3.1 ¿Podemos explicar que sucedió aquí?

3.4 Ejercicio

3.5 Bonus (Esto requiere investigar no basta con lo que aprendimos)

3.6 Recursos

4 RMarkdown

Placeholder

4.1 ¿Qué es R Markdown?

4.2 Librerías

4.3 Crear un R Markdown

4.4 Encabezado YAML

4.5 Sintáxis de Markdown

4.5.1 Lo esencial

4.5.2 Listas

4.5.3 Mostrar ejemplos de de código (sin ejecutarlos)

4.5.4 Figuras

4.6 MD -> HTML

4.7 Código R

4.8 Recursos

# load libraries
list.of.packages <- c("sf","mapview","terra","dplyr","RColorBrewer")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages,dep=TRUE,repos = "http://cran.us.r-project.org")
# library(landscapemetrics)
#library(raster)
library(RColorBrewer)
library(terra)
library(sf)
library(dplyr)

library(mapview)

5 Uso de R como Sistema de Información Geográfico

5.1 Clasificación de mapas y manejo de cartografía

Algunos insumos para manejo de cartografía que usaremos:

  • Mapas vectoriales de Chile Aquí
  • Procesamiento de datos geográficos con R Aquí o aquí

5.2 Usos de suelo para Chile

Importamos raster de cobertura de usos de suelo para Chile. Esta cobertura viene del lab de Jaime Hernandez en la U de Chile y esta documentada aqui y disponible aqui.

Luego de descargar los datos y descomprimirlos en nuestro directorio de trabajo podemos “cargar” el archivo TIF como un objeto en R

chiler <- terra::rast("/home/horo/Dropbox/tmp/coberturas/LC_CHILE_2014_b.tif")

Miremos la descripcion del archivo, fijemosnos en qué sistema de coordenadas está.

chiler
#> class       : SpatRaster 
#> dimensions  : 144511, 23799, 1  (nrow, ncol, nlyr)
#> resolution  : 30.00002, 29.99993  (x, y)
#> extent      : -8139.377, 705831, 3730098, 8065417  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) 
#> source      : LC_CHILE_2014_b.tif 
#> name        : LC_CHILE_2014_b 
#> min value   :               0 
#> max value   :            1210

Vemos que sus coordenadas están definidas en UTM 19S

Podemos también ver los datos con plot().

# ## sacar comentario para ver mapa
# plot(chiler,main="Usos de suelo")

5.3 Región de los Ríos

Esta es una cobertura en formato KML de la región de los ríos.

Fuente: aqui.

dir = 
reg = st_read("/home/horo/Dropbox/cursos/uach/CBIT241/practicos/Metricas2/Regional.shp") %>%
  filter(Region=='Región de Los Ríos')
#> Reading layer `Regional' from data source 
#>   `/home/horo/Dropbox/cursos/uach/CBIT241/practicos/Metricas2/Regional.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 17 features and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -12184470 ymin: -7554436 xmax: -7393642 ymax: -1978920
#> Projected CRS: WGS 84 / Pseudo-Mercator

plot(reg["Region"], axes = TRUE, graticule = TRUE)

5.4 Unificamos proyección geográfica

Ahora debemos reproyectar a UTM19 (epsg:32719)

r14_utm = reg |>
  filter(Region == 'Región de Los Ríos') |>
  st_transform(crs=32719)
plot(r14_utm["Region"], axes = FALSE, graticule = TRUE)

Podemos usar la provincia de Valdivia, del Ranco o de Chiloe en vez de toda la región

fuente: Provincias.zip

prov = read_sf("/home/horo/Dropbox/cursos/uach/CBIT241/practicos/Metricas2/Provincias.shp") %>%
  dplyr::select('Provincia') %>%
  dplyr::filter(Provincia == 'Chiloe') %>%
  st_transform(crs=32719)
#vald=st_transform(prov[prov$Provincia=='Valdivia',],crs=32719) # reproyeccion a UTM19S !!
#ranco=prov[prov$Provincia=='Ranco',]
par(mfrow=c(1,2))
plot(prov, axes = TRUE, graticule = TRUE)

#plot(ranco['Provincia'], axes = TRUE, graticule = TRUE)

Ojo que solo hemos reproyectado la provincia de Valdivia y no Ranco!

5.4.1 Comunas de Chile

#library(curl)
#com_tmp = tempfile()
#com = curl_download("https://www.bcn.cl/obtienearchivo?id=repositorio/10221/10396/2/Comunas.zip",com_tmp)
#unzip(com)

lr = read_sf("/home/horo/Dropbox/cursos/uach/CBIT241/practicos/Metricas2/comunas.shp") %>%
  dplyr::filter(Region == 'Región de Los Ríos'| Region == 'Región de Los Lagos') %>%
  dplyr::select(Comuna, Provincia, Region) %>%
  st_transform(crs=32719)

En R existen librerías que permiten la visualización interactiva como mapview

library(mapview)
mapview(lr,alpha=.3,zcol='Comuna')

5.4.2 Clipping y reclassificación de usos

Debemos (re)cortar el raster con el poligono reproyectado para trabajar unicamente con área de interes. Para eso, usamos funciones de la librería terra. La función crop recorta a la extensión de chiler, pero solamente puede recortar coberturas de tipo raster en formas rectangulares, pues estamos tratando con una matriz (!). Para eliminar todo lo que no está en nuestro polígono de referencia usamos mask, lo que va a enmascarar lo que no está en nuestro polígono de referencia. Finalmente, trim elimina columnas, o filas, que solo tengan NaN.

La idea es tener una cobertura raster lo mas pequeña posible para conservar memoria. Por lo mismo, se eliminan las coberturas/objetos que no usaremos mas con la función rm.

la_comuna = 'Máfil'
v_paisa <- subset(lr,Comuna == la_comuna) # el paisaje que vamos a estudiar (en vectorial)


#LC <- chiler
LC <- terra::crop(chiler,vect(v_paisa) )
LC <- terra::mask(LC, vect(v_paisa) )
LC <- terra::trim(LC)
#plot(LC, axes = FALSE, background="grey",mar=c(2,2,2,6))
plot(LC, axes = FALSE,mar=c(2,2,2,6))


#rm(chiler,r14,r14_utm)

5.4.3 Reclasificamos a usos que nos interesan

Vamos a ver qué clases de usos existen. Ver equivalencia en documentación. Po ejemplo, los valores entre 200 y 300 son alguna categoría de bosque (ver mas abajo).

sort(unique(values(LC)))
#>  [1]  120  130  140  150  211  212  241  251  311  312  320
#> [12]  330  410  420  510  610  630  640  800  920  931  932
#> [23] 1210

Es mas informativo ver la frecuencia de cada uso con un gráfico de barras


barplot(LC,main='Número de pixeles por clase',horiz=TRUE,las=1)

Ahora necesitamos reclasificar los Usos de Suelos (Land Covers). Usaremos la Tabla 1 de la documentación.

5.4.4 Nos interesan retener los siguientes usos:

  • Oceano : 0
  • Cultivos: 1 a 200
  • Bosque Nativo: 210,220 y 230)
  • Plantaciones (240 y 250)
  • Pastizales (300)
  • Matorrales (400)
  • Humedales (500)
  • Cuerpos de agua (600)
  • Otros (800,900,1000 y 1200)

Para reclasificar, generamos una tabla con los valores de reclasificacion

#clases que nos interesan
clases_usos = c("Cultivos","Bosque Nativo","Plantaciones","Pastizales","Matorrales","Humedales","Agua","Otros")
# matriz para reclasificar, la primea columna es desde, la 2da hasta y la 3ra el valor de clase que le daremos (ver: ?terra::classify)
reclass = c(
  -Inf,1,NaN,
  1,200,0,  # Cultivos
  200,240,1, # Bosque Nativo
  240,300,2, # Plantacion
  300,400,4, # Pastizales
  400,500,5,
  500,600,6,
  600,800,7,
  800,2000,8
)
reclass_m = matrix(reclass,ncol=3,byrow=TRUE)

col2 = c("red","darkgreen","brown","yellow","wheat","lightblue","blue","white")

LC2 = terra::classify(LC,reclass_m)
usos = data.frame(id=sort(unique(values(LC2))), usos=clases_usos)
levels(LC2) = usos # establece los valores del raster como categórico, como clases

par(mar=c(2, 5, 2, 1))
barplot(LC2,main='Número de pixeles por clase',names=clases_usos,horiz=TRUE,las=1,col=col2,cex.names=.6)

#rm(LC,reclass_m, reg)

Luego, vamos a visualizar el raster clasificado

plot(LC2, type="classes",col=col2,main=la_comuna, axes=F, mar=c(1,1,2,6),legend=T)

Podemos ver también solo aquellos parches de bosque nativo

par(mfrow=c(2,2))
plot(LC2 ==1,main="bosque nativo", axes=F,legend=F,mar=c(1,1,3,1),col=c("grey","darkgreen"))
plot(LC2 ==2,main="plantación", axes=F,legend=F,mar=c(1,1,3,1),col=c("grey","brown"))
plot(LC2 ==4,main="pastizales", axes=F,legend=F,mar=c(1,1,3,1),col=c("grey","yellow"))
plot(LC2 ==0,main="cultivos", axes=F,legend=F,mar=c(1,1,3,1),col=c("grey","red"))